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ABSTRACT 

We describe a simple, efficient, robust and fully automatic algorithm for the determination of 
a Multi-Gaussian Expansion (MGE) fit to galaxy images, to be used as a parametrization for 
the galaxy stellar surface brightness. In most cases the least-squares solution found by this 
method essentially corresponds to the minimax, constant relative error, MGE approximation 
of the galaxy surface brightness, with the chosen number of Gaussians. The algorithm is 
well suited to be used with multiple resolution images (e.g., Hubble Space Telescope [HST] 
and ground-based). It works orders of magnitude faster and is more accurate than currently 
available methods. An alternative, more computing intensive, fully linear algorithm, that is 
guaranteed to converge to the smallest solution, is also discussed. Examples of MGE 
fits are presented for objects with HST or ground-based photometry, including galaxies with 
significant isophote twist. 

Key words: galaxies: kinematics and dynamics - galaxies: structure - techniques: image 
processing - stellar dynamics 



1 INTRODUCTION 

Digital images of galaxies which are small by today's standards, 
may still contain over one million pixels. It is natural to try to distill 
the information contained in this large number of pixels into a set 
of quantities that can be used easily and efficiently to derive infor- 
mation on the physical properties of the observed object. The most 
common approach currently consists in fitting ellipses of increasing 
semimajor axis to the galaxy images (e.g., Carter 1978; Kent 1984; 
Lauer 1985; Jedrzejewski 1987; Franx, Illingworth & Heckman 
1989; Peletier et al. 1990). Intensity profiles, ellipses shape and 
deviation from ellipses, parametrized by a few Fourier terms, are 
obtained with these methods. A nice feature of this parametrization 
is that it gives a physically meaningful description of the galaxy 
photometry in terms of the ellipticity, diskyness and boxyness at 
each radius. 

One limitation of the ellipse fitting methods however is that 
strong deviations of the isophotes from ellipses cannot be easily 
modeled with a small number of Fourier terms. This makes these 
methods not well suited to the photometric modeling of multicom- 
ponent objects such as lenticulars and spirals, which have a bulge, 
a main disk and often an embedded nuclear disk or a bar. More 
importantly, if the deviations from ellipses are non negligible it be- 
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comes non trivial to use the fitted isophotes for the construction of 
realistic galaxy dynamical models due to the complexity of depro- 
jection and the evaluation of the gravitational potential. 

In this paper we focus on a completely different approach to 
the photometric modeling, the Multi-Gaussian Expansion (MGE) 
method, that overcomes both of these limitations. The original 
idea of the application of the MGE method to the problem of 
the deconvolution and deprojection of galaxy images comes from 
Bendinelli (1991). This idea was generalized to the non-spherical 
case and made applicable to real galaxies by Monnet, Bacon & 
Emsellem (1992) and further developed by Emsellem, Monnet & 
Bacon (1994a) and Emsellem, Dejonge & Bacon (1999). 

The MGE method consists of a series expansion of galaxy im- 
ages using two-dimensional (2D) Gaussian functions. The Gaus- 
sians have the beneficial properties that both the convolution (e.g., 
to take seeing or PSF effects into account) and the deprojection 
(to derive the intrinsic stellar luminosity density from the observed 
galaxy photometry) can be performed analytically in a simple and 
efficient manner. As shown by Emsellem et al. (1994a), many other 
dynamical and photometric quantities can be evaluated easily and 
accurately when the density is expressed in MGE form. For exam- 
ple, the MGE potential can be computed with a single integration, 
as opposed to the two that are required when the intrinsic density is 
stratified on similar triaxial ellipsoids, and three in the general case. 
In the case of ,f{E, L^) axisymmetric MGE dynamical models the 
velocity moments predicted from the leans equations, already pro- 
jected onto the sky, can be expressed with only a double integra- 
tion. The even part of the distribution function can also be easily 
retrieved from an MGE density distribution via the Hunter & Qian 
(1993) formaUsm. 
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The MGE parametrization is one of the few simple parame- 
trizations that are general enough to reproduce the surface bright- 
ness of realistic multicomponent objects (e.g., spirals with multiple 
disks). It has already been used for the modeling of a number of 
galaxies (Emsellem et al. 1994b, 1996, 1999; Emsellem 1995; van 
den Bosch et al. 1998; van den Bosch & Emsellem 1998; Cretton 
& van den Bosch 1999). 

What is however still missing from the MGE machinery is an 
accurate, easy to use, generally available, robust and automatic al- 
gorithm for the determination of an MGE expansion from galaxy 
images. In this paper we present such an algorithm. The method we 
describe here brings the MGE fitting phase of a dynamical model- 
ing process to the same level of the LOSVD extraction phase (e.g., 
van der Marel & Franx 1993) in terms of robustness and ease of 
use. 

This paper is organized as follows. In Section 2 we define our 
notation and give some useful formulae for the application of MGE 
models to the dynamical modeling of galaxies. In Section 3 we dis- 
cuss our new MGE fitting method. In Section 4 we give information 
on the availability of the software implementing the methods dis- 
cussed in this paper. Some conclusions are given in Section 5. 



2 MGE FOIfMALISM 

In this section we rederive, in a slightly different notation, results 
from Monnet et al. (1992) and Emsellem et al. (1994a) that are usu- 
ally needed for the construction of dynamical models, starting from 
the MGE photometric models we describe in this paper. The case 
of Gaussians with the same centre and circular PSF is considered 
here, while the general case is discussed in detail in Emsellem et 
al. (1994a). 

Let {x',y',z') be a system of coordinates centred on the 
galaxy nucleus, with the z' axis pointing towards the observer. The 
MGE projected surface brightness can be written as 



E{R',e') 



N 



■ exp 



2-1 



y/2 



with 



x'j = R'sin{e' -ipj), 
j/^ =i?'cos(e'-V,), 



(1) 



(2) 



where {R',9') are the polar coordinates on the plane of the sky 
{x',y'). Here A'^ is the number of the adopted Gaussian compo- 
nents, having total luminosity Lj, observed axial ratio < < 1, 
dispersion aj along the major axis, and position angle (PA) ^j, 
measured counterclockwise from the y' axis to the major axis of 
the Gaussian. 

Given that the convolution, deprojection and potential calcu- 
lations, can be carried out separately for each single Gaussian com- 
ponent of equation ([l]), in what follows the index of the Gaussian 
is not written explicitly and only one Gaussian is considered. The 
results for the complete MGE model are obtained by summing the 
contributions over the A'^ Gaussians components. 

2.1 PSF convolution 

The MGE surface brightness has to be convolved with the instru- 
mental PSF before comparison with the observed surface bright- 
ness. Assume the normalized PSF can be written as the sum of M 
circular Gaussians 



PSF(i?') = 



2'!vaf 
fc=i " 



exp 



(3) 



with X/fc=i ~ 1- '^"'^''1 luminosity L of one Gaussian, 
having surface brightness E, does not change after convolution 
E = E ® PSF with the PSF. The convolved Gaussian can then 
be written explicitly as 



E{R',e') = L 
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■ exp 
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2.2 Deprojection 

The deprojection of a galaxy surface brightness distribution is gen- 
erally non unique (e.g., Rybicki 1986; Franx 1988). The MGE 
deprojection discussed here represents a possible solution to the 
problem. It generally produces smooth and natural-looking in- 
trinsic densities, consistent with the observed photometry. The 
nonuniquness of the triaxial deprojection on a set of real galaxies 
and its effects on the dynamical modeling will be the subject of a 
coming paper (Cappellari & de Zeeuw, in preparation). 

2.2.1 Triaxial case 

If a good MGE fit to the photometry can only be obtained by al- 
lowing the PA tpj of the Gaussians to be different from each other, 
then the object can not be axisymmetric. The following deprojec- 
tion method can be used to produce an intrinsic density model con- 
sistent with the observations. 

Consider one single Gaussian component taken from equa- 
tion (|l|). Since the isophotes of this component are similar ellipses 
it can be deprojected by assuming the intrinsic density is strati- 
fied on triaxial ellipsoids (Contopoulos 1956; Stark 1977; Binney 
1985). The intrinsic luminosity density of this component can be 
written again as a Gaussian with the same total luminosity L 



p(x,y,z) 



■ exp 



1 

2^ 



y 



+ 



(6) 



(\/27r (T)^pg 

Here (x, y, z) is a system of coordinates as in Binney (1985), cen- 
tred on the Gaussian centre and aligned with its principal axes. The 
axes are assumed to be such that the longest axis of the Gaussian 
lies along the x-axis and the shortest axis along the z-axis, so that 
< g < P < 1- The two usual polar coordinates {6, 0) define the 
orientation of the line-of-sight with respect to the principal axes of 
the object. Another angle il) is required to specify the rotation of 
the object around the line-of-sight, and uniquely define the orien- 
tation in space of the [x, y, z) coordinate system. This angle is 
defined by the requirement that the z axis project onto the y' axis, 
or equivalently that the x' axis lies in the {x, y) plane. 

Once the viewing angles {9, 4>, '4') have been assumed, given 
the observed axial ratio q' for the Gaussian, the intrinsic axial ratios 
p and q can be found using relations valid for general ellipsoidal 
bodies (e.g., de Zeeuw & Franx 1989) 

5' [2 cos 2i/;+sin 2^(sec 9 cot — cos 6 tan 0)] 



1-g 



2 sin^ 9 [5' cos ip{cos^-\-col <f) sec 9 sin — 1] ' 
5' [2 cos 2'!/)-|-sin 27/;(cos 6 cot 4>— sec 9 tan </>)] 
2 sin^ 9 [5' cos -i/j (cos 7/)+ cot 4> sec 9 sin 1/') — 1] ' 



(7) 



(8) 
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where 5' = 1 — q'^. It should be noted that in general a solution 
may not be found for all assumed viewing angles (e.g., Bertola et 
al. 1991). If a triaxial deprojection is sought for the galaxy, then the 
viewing angles have to be the same for all the Gaussians and the 
permissible viewing angles for the object, assumed to be triaxial, 
will lie in the intersection of the permissible viewing angles of each 
Gaussian (cf. Williams 1981). 



2.2.2 Axisymmetric case 

If an acceptable MGE model can be obtained by fixing the PA tpj of 
all the Gaussians to the same value tp, then the object can be depro- 
jected by assuming it is axisymmetric (triaxial solutions however 
can not be excluded [cf. Franx 1988]). In this case the formulae of 
the previous section simplify considerably. 

In the oblate axisymmetric case (p = 1), once an inclination 
i > (for i — the deprojection is degenerate) has been assumed 
for the galaxy (i = 90° corresponding to edge-on), one has: 



where 5 — l—p^,e = 1 — g^,Gis the gravitational constant and T 
is the mass-to-light-ratio. For oblate axisymmetric models (p = 1) 
this reduces to 



■dr. (12) 



In the spherical case (p = q — l)lo 



GTL 
R 



erf 



(13) 



(9) 



Already for 8a ^ R the error function erf [_R/(V2 a)] = 1 to 16 
digits accuracy and the spherical MGE potential can be replaced 
with high precision by the potential of a point mass enclosing the 
total mass A4 = TL of the Gaussian component. 

Since the a of the Gaussians in an MGE fit span various or- 
ders of magnitude in size, asymptotic expressions can also be con- 
veniently used to estimate the triaxial or axisymmetric potential of 
equations ( |ll| ) and ( |l^ at very small and very large radii as done 
in Cappellari et al. (2002). 



while in the prolate axisymmetric case (p — q) 



2 

Q = 



COS^ I 



(10) 



The oblate MGE deprojection above is only defined if 
cos^ i < q'j^ for all Gaussians. This means that the flattest Gaussian 
in an MGE fit dictates the minimum possible inclination for which 
the MGE model can be used in a dynamical model. However since 
the Gaussian components don't necessarily have physical signifi- 
cance one don't want to base any conclusion of a dynamical model 
on this minimum inclination. 

One simple way to avoid this problem consists in trying to 
increase the minimum axial ratio q' allowed in the MGE fit until 
the increases significantly or the galaxy contour cannot be ac- 
curately reproduced any more with any number of Gaussians. This 
(Jmax will provide a good estimate of the minimum inclination for 
which the given galaxy, assumed to be axisymmetric, can be depro- 
jected. 

Independently from an actual MGE fit, the limit to the above 
9max is set by the axial ratio q' of the ellipse having the same cur- 
vature of the galaxy isophotes along the major axis, at the point of 
maximum curvature. In fact there has to be at least one Gaussian 
flat enough to reproduce the maximum curvature along the major 
axis. 



2.3 Gravitational potential 

The potential generated by a triaxial Gaussian mass distribution of 
the form given in equation (^, can be evaluated using the classical 
Chandrasekhar (1969) formulae for densities stratified on similar 
concentric ellipsoids (see Binney & Tremaine 1987, p. 61). In the 
MGE case the potential can be written as a one-dimensional (ID) 
integral 



exp 



^(1-5T2)(1-.T2) 



dT, 



(11) 



3 MGE FITTING 

3.1 The brute-force 2D approach 

We start by considering the MGE fit to galaxies without isophote 
twists. The determination of an axisymmetric MGE parametriza- 
tion from the photometry, namely the fit of a sum of Gaussians with 
the same centre and PA, to a galaxy image seems to be a very sim- 
ple problem that can be efficiently solved with general nonlinear 
minimization algorithms. 

We assume for simplicity that the PA i/i of all Gaussians is 
known. We don't fit ^ since it is very easily measured before fitting. 
We found it fast and accurate to determine it from the weighted sec- 
ond moments of the surface brightness above a given level, as done 
by automatic galaxy extraction packages (e.g., Bertin & Amouts 
1996). It may also be measured with standard photometric pro- 
grams, or may be dictated by other requirements in the dynamical 
model (in general a more accurate is required in the centre where 
the kinematical data are available). We will also always assume that 
the MGE models are analytically convolved with a circular MGE 
PSF before comparing with the data using Equations and (^. 

The straightforward MGE fit approach consists of fitting the 
nonlinear function of equation (|l]) to the pixels of the image by 
minimizing the (pixel integration can be taken into account in 
the very innermost few pixels) 



E 



~ ^(s^i) Vj] 



(14) 



where C^j are the pixel counts, Ai.j is a corresponding weight 
factor, depending on the photometric and systematic errors, and the 
sum is done over all (good) image pixels. Positivity constraints are 
generally enforced on the Lj, although in principle only the total 
surface density needs to be positive everywhere. 

Although the convergence can be improved by starting the fit 
with only a few Gaussians and then adding new components as long 
as the residuals decrease appreciably (Emsellem et al. 1994a), we 
have discovered that this brute force approach does not work well 
in practice for two main reasons: 
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(1) the nonlinear function above iias a large number of lo- 
cal stationary points. In fact suppose that we have found a mini- 
mum of the with A'' Gaussians. A new configuration with A'^ + 1 
Gaussians, where the new Gaussian is equal to one of those previ- 
ously found, is a saddle point for the . Consequently even state- 
of-the-art nonlinear fitting algorithms have a high chance of stop- 
ping or slowing down considerably in a narrow valley in the SA'^- 
dimensional space (with ~ 10 — 15) close to one of these many 
stationary points, unless the initial guess is very near the true mini- 
mum; 

(2) the "correct" weighting of the pixel measurements tends to 
produce a model that underestimates the nuclear density due to the 
fact that the total of equation ( |l4| ) is dominated by the contribu- 
tion of many pixels at large radii. This is not what is needed for an 
accurate dynamical model of the nucleus. 

The net result is that in the above approach the fit proceeds by 
trial and error, with a lot of manual intervention. An automatic fit 
generally takes many hours on a 1 GHz PC and usually converges 
to a solution that can still be improved manually. 



3.2 An efficient ID algorittim 

To solve the convergence problems of an MGE fit we start by 
analysing the peculiarities of a fit to a ID profile, as can be obtained 
by measuring the photometry of a galaxy with constant elliptic- 
ity. Since galaxy profiles are generally well described by multiple 
power-law regimes (e.g., Lauer et al. 1995), in this section we start 
by trying to satisfy the natural requirement that the fit has to "look 
good" on a log-log plot. To obtain this result we do the following: 

(1) we sample the profile logarithmically in radius. This is usu- 
ally done in galaxy photometry to improve the S/N at large radii 
while still maintaining all the information at the smaller radial 
scales; 

(2) we require the fit to optimise the relative error, namely we 
minimize 



2 



E 



(15) 



where M is the number of photometric data points Cj in the profile, 
Rj is the corresponding isophote major axis and E(_R) is the ID 
version of equation ([l]). 

The reason why Gaussians are so successful in reproducing 
power-law profiles lies in the fact that every Gaussian provides a 
very localized contribution to the profile, only for radii close to 
its (Tj. In fact I) at small radii {R <^ aj) Gaussians become very 
rapidly constant, while power-laws continue to rise, 2) at large radii 
(R ^ (jj ) Gaussians go to zero much more rapidly than a power- 
law. It appears clear that the MGE best-fitting solution for a loga- 
rithmically sampled power-law profile consists of a series of Gaus- 
sians that are also logarithmically spaced, both in aj and in height, 
within the ranges imposed by the data (see Fig.[l|). 

We now show that this solution is the true global minimum. In 
fact our least-squares MGE solution, with positive Gaussians, for a 
given decreasing function of radius, in a given interval, is charac- 
terized by the condition that the maximum relative error is reached 
precisely 2N + 1 times, with alternating signs (bottom panel of 
Fig. |l]). A similar condition also exists for the unique minima:^ 



The solution minimizing the maximum absolute error. 



10 



-4 



10 



0.01 





0,01 0.1 1 '0 100 

Figure 1. ID MGE fit to the power-law profile . Top panel: the pro- 
file and the N = Id best-fitting Gaussian components ai'e shown. Bottom 
panel: radial variation of the relative error in the profile fit. Note that the 
maximum error is reached precisely 2N + 1 times, with alternating signs. 
The vertical lines indicate the Cj of the Gaussians shown in the upper panel. 



solution of other approximating functions e.g., Chebyshev polyno- 
mials and rational functions (see e.g.. Press et al. 1992). The fact 
that the minimax condition is satisfied for our MGE fits shows that 
our algorithm is able to converge to the minimax MGE approxi- 
mation of a power-law. The same minimax condition can be veri- 
fied in the case of MGE fits to more realistic profiles composed of 
multiple power-law regimes (Fig. and essentially for any func- 
tion that appears smooth in a log-log plot. Here the Gaussians will 
not be precisely logarithmically spaced in the Oj , but will be more 
closely spaced where the profile is steeper (Fig. Also in these 
cases however the least-squares solution essentially provides the 
minimax MGE approximation of the given function. 

We emphasize the fact that the logarithmic sampling of the 
profiles is a necessary condition for the least-squares solution to 
coincide, within the numerical approximations, with the minimax 
solution for the power-law. A linear sampling would have produced 
an MGE fit with smaller errors at large radii and larger errors to- 
wards the centre. 

With real photometric data the above condition for the mini- 
max solution cannot be always verified, due to the noise and to the 
fact that the observed profile is known at a finite number of points. 
The convergence to a global solution however can still be recog- 
nized, by simple visual inspection of the result, from the fact that 
every Gaussian has to give a significant contribution to a specific 
part of the profile: a Gaussian that remains always lower than any 
other Gaussian is generally an indication that the solution is not the 
global minimum (unless the profile becomes very flat). 

To determine the fits of Fig. we made use of the important 
fact that equation (|l|) is linear in the Gaussian total luminosity Lj . 
Since a global solution for the linear variables Lj can be found 
easily for any given combination of the nonlinear ones (jj, it is 
crucial to separate the treatment of these two sets of variables. In 
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Figure 2. Same as in Fig. |l]for the triple power-law profile of the intrinsic 
density of M32, as determined by van der Marel et al. (1998). With N = 16 
Gaussians the maximum error is 0.7 per cent (excluding border effects) in 
the range (X'Ol-SOO", over ten orders of magnitude in the density. 
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Figure 3. Same as in Fig. |l]for the double power-law profile with flat core 
(1 + R)^'^. The maximum error is 0.3 per cent (excluding border efi'ects) 
over ten orders of magnitude. 

practice our ID MGE fitting algorithm consists of the following 
steps: 

(i) we start with the CTj logarithmically spaced in radius between 
the maximum radius of the measured photometric data J?max and a 
minimum radius J?min = 0.75 o"psf. We know this is a good guess 
for general multiple power-law profiles; 



(ii) we solve the nonlinear least-squares minimization prob- 
lem for the logfjj, imposing the constraints that 7?min < o"j < 
i?max. The actual computation is done using the robust Levenberg- 
Marquardt implementation due to More, Garbow & Hillstrom 

(1980)|[J| The adopted implementation accepts bound constraints on 
the variables and is able to deal with degenerate Jacobians. The lat- 
ter is a necessary feature since some of the variables Lj may be 
zero in the final solution; 

(iii) for every choice of the aj variables, during the above itera- 
tive optimisation process, we find the best-fitting value for the Lj, 
with the constraints Lj > 0. This is a Non-Negative Least-Squares 
(NNLS) problem that can be efficiently solved using the Lawson & 
Hanson (1974) algorithm^. 

In case a solution with negative Gaussians is needed (e.g., to 
reproduce a profile that is not monotonically decreasing with ra- 
dius) the solution of the linear system with NNLS, in step (iii) 
above, must be replaced by a solution using the Singular Value De- 
composition (SVD) to provide a zero order (i.e. minimum norm) 
regularization on the solution (e.g.. Press et al. 1992). Other linear 
regularization constraints could also be enforced on the solution. 
A direct non-regularized solution of the linear system is extremely 
ill-conditioned and will produce a result containing elements Lj of 
very large absolute value and different sign. An MGE model ob- 
tained in this way presents the risk of cancellation of significant 
digits when deriving any physical quantity. 

Fig. ^ shows that the maximum relative error in an MGE 
fit to multiple power-law profiles decreases almost exponentially 
with increasing number of Gaussians, approximately as fast as 
e oc e^" *^. This is not surprising, given that every Gaussian di- 
verges exponentially from the power-law starting from the point of 
closest approach. Although the actual numbers will depend on the 
profiles and on the adopted radial ranges, these examples provide 
an estimate of the number of Gaussians that are required to reach a 
given accuracy in the fit, and explain why 10-15 Gaussians gener- 
ally produce a very good MGE fit to realistic Galaxy profiles. 

We finally mention that an asymptotic power-law nuclear cusp 
can also be imposed on the MGE profile as described in Emsellem 
et al. (1999). Although this choice breaks the elegance of the MGE 
parametrization somewhat, this detail does not change in any way 
the essence of our fitting algorithm. 



3.3 Extension of the ID algorithm to 2D 
isophote twist 



images without 



We now extend the ID algorithm described above to the fit of actual 
2D images. We perform this extension by starting from the obvious 
requirement that the 2D algorithm provides exactly the same results 
of the ID case (to be precise the same Lj and Oj) when it is used to 
fit a 2D image that presents the same major axis profile of the ID 
case and has perfectly elliptical isophotes with constant axial ratio 
q and PA. 

This goal is achieved by fitting "in parallel" the MGE model 
of equation (Q), in polar coordinates, to a certain number A'^bcc of 
photometric profiles, measured along sectors uniformly spaced in 



We have used an IDL por ting of the MINPACK-1 code bv Craig B. Mark- 



wardt that is available from tittpV/cow.physics. wisc.edu/~craigm/idl 

We have ported to IDL their 1995 Fortran 90 version of the Bounded- 



Va riables Least-Squares (BVLS) code, a generalization of NNLS, available 
on ^ittp://www.netlib.org 
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Figure 4. Variation of the maximum relative error in the MGE fit of Fig. ^ 
(solid line) and in the fit of Fig. |^ (dashed line) as a function of the number 
of Gaussians adopted in the expansion. 



angle from the major axis to the minor axis (sectors in the four 
quadrants are averaged together before fitting). All profiles include 
the central pixel and proper integration of the Gaussians on that 
pixel was performed for an accurate data-model comparison. The 
sampling of the photometry along the sectors is now logarithmi- 
cally spaced in the elliptical radius m'^ — x'^ + y'"^ /q'^ , where q' 
is a representative axial ratio of the galaxy isophotes (see Fig. ^. 
This latter detail ensures that the requirement stated in the previous 
paragraph is precisely satisfied also for flattened objects. This way 
of sampling galaxy photometry is a time proven practice common 
to most photometry packages (e.g., Jedrzejewski 1987). In the fit 
the same relative error is assigned to all measurements as in the ID 
case, following equation (|l5|). 

The actual photometric measurements are obtained by select- 
ing sets of pixels on the image, based on their angular position and 
logarithmic radius ranges. Once the pixels have been selected, the 
luminosity weighted average coordinates of each photometric point 
are evaluated from the pixels themselves and not from the original 
sectors. This ensures that the coordinates of sectors at the border, 
or containing masked pixels, are properly computed and, even more 
importantly, it permits an accurate modeling of the galaxy surface 
brightness distribution close to the centre, where the discrete nature 
of pixels cannot be neglected. In particular no surface brightness in- 
terpolation is performed in the centre, but the observed pixel values 
are used in the fit. To guarantee a proper sampling of the image, 
both in the nucleus and at large radii, 24 sectors have been adopted 
for every decade in radius in pixels units (this corresponds to the 
usual factor 1.1 separation between successive radial intervals). 

With these choices, when fitting a galaxy with slowly vary- 
ing ellipticity with radius (as is often the case for single component 
objects like ellipticals), if the input guess q' for the axial ratio is 
close to the galaxy average value, the procedure essentially reduces 
to A'sec ID similar fits and converges very quickly by varying al- 
most only the (jj as done in Section 3.2. If the axial ratio is not 



constant or the isophotes are not ellipses (as in the case of spirals), 
both the Oj and q'^ will have to be optimised together. In practice 
we still apply steps (i)-(iii) of Section 5.2, with the only difference 
that we also fit the nonlinear parameters q'j in step (ii). To make 
the fit more robust, bound constraints are given on the variables: 
Rmin < Oj < i?max and Q < q'j < 1. Experience on a number of 
complex objects has shown that the above procedure will still con- 




Figure 5. This image illustrates the sectors that have been used to produce 
an MGE model for the WFPC2/F814W image of the flattened SO galaxy 
NGC 4342. In the fit A^scc = 19 sectors have been used, spaced in angle 
by 5° and covering the whole galaxy from the major to the minor axis. In 
this figure only one every thi'ee sectors is shown. These ai'e the sectors from 
which the profiles shown in Fig. ^ have been extracted. Note that the radial 
sampling along the sectors follows a characteristic isophote of the galaxy. 
The 34" X 34" PC 1 image of the galaxy is visible in the background. 



verge very quickly and automatically to an excellent MGE fit to the 
photometry, since the initial guess for the solution is already good 
enough to permit the Levenberg-Marquardt optimisation algorithm 
to converge to the global minimum. 

As an example we present in Figures I the profilesj^] of an 
MGE fit to the elliptical galaxy M32 (cf. Bendinelli et al. 1992). 
The model was fit simultaneously to a WFPC2/F814W image and 
an /-band ground-based one (Peletier 1993). The A'^ = 11 Gaus- 
sians MGE fit for this galaxy was obtained with Ascc = 19 sectors, 
spaced by 5°, covering the whole galaxy, from the major to the 
minor axis. As expected the Gaussians succession is very regular 
and predictable, as in the ID MGE fit previously shown in Fig. ^. 
The comparison between the contours of the model and the data is 
presented in Fig. ^ Similarly to what was done by Emsellem et al. 
(1994a), to determine the number of Gaussian components in the 
fit we start e.g., from A'^ = 10 Gaussians and keep adding Gaus- 
sians and restarting the fit, until the best-fitting stops decreasing 
appreciably (e.g., less than 1 per cent). 

Instead of convolving the MGE model with the different PSFs 
during the fit we simply exclude the radial range that is affected 
by PSF effects, of the lower resolution image, and only convolve 
the model with the HST PSF. This is done to avoid spoiling the 
fit to the high resolution HST image by any PSF mismatch in the 
ground-based data. Sky subtraction and relative flux calibration of 
the frames is standard for the MGE method as for any other photo- 
metric modeling (e.g., Franx et al. 1989) and will not be discussed 



tt A value Rmin = [Vs + sinh~^(l)]/6 ~ 0.38 pixels has been as- 
signed in the plots to the radius of the central pixel, corresponding to the 
unweighted average radius measured from the centre of the pixel. 
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Figure 6. Le^ Panels: Comparison between the WFPC2/F8 14W + ground- 
based 7-band photometry of M32 (open squares) and the corresponding 
A'^ = 11 Gaussians MGE best-fitting model (solid line). The individual 
convolved Gaussians components are also shown. From top to bottom the 
profiles are measured along 5° wide sectors, linearly spaced in angle be- 
tween the major (0°) and the minor axis (90°). Only 7 representative sec- 
tors of the of the A^soc = 19 sectors actually used in the fit are shown here. 
Right Panels: radial variation of the relative error along the profiles. The 
RMS error is 2.2 per cent. Most of this scatter is due to the resolution of 
M32 into stars. 



here in detail. The sky subtraction for the ground-based image was 
determined by requiring the outer profile to fall of as a power-law; 
the sky level and relative calibration of the HST image was adjusted 
to ensure agreement in the overlap region. 

An MGE fit to the more "difficult" SO galaxy NGC 4342, 
which has a small nuclear stellar disk, in addition to a main outer 
disk (van den Bosch et al. 1998) is shown in Fig. |8| Also in this 
case the MGE fit with N = 13 Gaussians to the WFPC2/F814W 
image of this galaxy was obtained with A'scc = 19 sectors, spaced 
by 5°, covering the whole galaxy, from the major to the minor axis. 
A comparison between the contours of our MGE model and of the 
observed image is shown in Fig. ^ 

In Fig. |l^ is presented an example of an MGE fit to the Sa 
galaxy NGC 4698, which has non-convex isophotes (Bertola et al. 
1999; Cappellari 2000) and for this reason necessarily cannot be 
reproduced by a sum of positive Gaussians. In this case the posi- 
tivity constraint for the Lj has been eliminated and the linear part 
of the Gaussian fit has been performed using SVD as described in 
Section 3.2. 




• (/ 



Figure 7. Top panel: Contour map of the central 34" X 34" of the (1280 s) 
WFPC2/F814W (/-band) image of M32. The noisy appearance of this high 
signal-to-noise HST image is due to the resolution of the galaxy into stars. 
Bottom panel: Same as in upper panel for the 200" X 300" ground-based 
(1 s) INT 7-band image. Superimposed on the two plots are the contours 
of the intrinsic MGE surface brightness, convolved with the proper PSF. In 
this plot, as in the following ones the contours are logarithmically spaced, 
but otherwise arbitrary. 
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Figure 8. Same as in Fig. |^ for the N = 13 Gaussians MGE fit to tlie 
WFPC2/F814W pliotometry of NGC 4342. Tlie RMS error is 1.8 per cent. 
Part of the data-model discrepancy at large radii along the major axis is due 
to a non perfect flatfielding of the PC 1 CCD close to the border. 



NGC 4698 has a disk with prominent dust lanes. The tech- 
niques that can be used to determine an MGE model of a dusty 
galaxy are almost the same one can use to correct standard pho- 
tometric measurements for dust effects (e.g., Carollo et al. 1997). 
While performing an MGE fit one can mask (exclude from the fit) 
the galaxy regions that are affected by clumpy dust absorptions. 
This method clearly works better if some symmetries can be as- 
sumed for the underlying galaxy. This was done to obtain the MGE 
model of Fig. |l^ with constant PA. MGE models like this can also 
be used to estimate and correct for dust absorption effects (Em- 
sellem 1995). When images in more than one band are available, 
a colour excess map [e.g., E{V — I)] can be used to perform a 
first-order correction for dust absorption as done by Cappellari et 
al. (2002) to determine an MGE model of IC 1459. 

The MGE fitting method discussed in this paper was already 
successfully applied to another sample of ten galaxies in Cappel- 
lari (2000). The method was also tested on twenty elliptical and 
lenticular galaxies with HSTAVFPC2 photometry, taken from the 
SAURON representative sample (de Zeeuw et al. 2002). These 
models will be presented elsewhere. We have generally found that 
the RMS error of our MGE models is dominated by the noise or by 
intrinsic asymmetries in the photometry, and is of the order of ~ 1 
per cent or less. In most cases a good fit is obtained in about one 




Figure 9. Contour maps of the (200 s) WFPC2/F814W (/-band) image of 
NGC 4342 at two different scales: 8" X 8" (top panel) and 34" X 34" 
(bottom panel). Superimposed on the two plots are the contours of the in- 
trinsic MGE surface brightness (Fig. convolved with the WFPC2 PSF. 
This figure can be compared with Fig. 1 in Cretton & van den Bosch (1999). 



minute on the same 1 GHz machine that took many hours to obtain 
a poorer fit with the brute force approach. 

We notice that the adopted MGE fit approach also has the ad- 
vantage that it is very simple to fit together multiple photometric 
profiles obtained from different images. It is also easy to combine 
high resolution images with published photometric profiles at large 
radii. In fact during the fit the program only needs to know the 
measured values of the surface brightness and the corresponding 
polar coordinates, and does not have to deal with multiple images 
directly. 
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Figure 10. Contour plot of a V-band ground-based image of the Sa 
galaxy NGC 4698, taken with the Vatican Advanced Technology Telescope 
(VATT). Overlaid are the contours of the best-fitting MGE model for this 
galaxy, obtained without imposing the usual positivity constraints on the 
Gaussians. Only in this way also the non-convex central isophotes can be 
reproduced. The prominent dust lanes in the bottom part of the image and 
the bright stars have been masked during the MGE fit. 



3.4 A fully linear alternative 2D algorithm 

An alternative approach to the 2D MGE fit of an image, using Gaus- 
sians with the same PA and the same centre, consists of transform- 
ing the whole problem into a NNLS problem, that can be solved 
with the ad-hoc Lawson & Hanson (1974) algorithm. This can be 
done by generating a large grid of all possible Gaussians that may 
end up in the final solution. We sample logarithmically with 
Rmin < Oj < J?max; for cvery choice of Gj we sample q'^ linearly 
with < < 1. The non-negative superposition of the Gaussians 
that best fit the data is then found with NNLS. This approach is 
memory intensive since a copy of the data has to be kept in mem- 
ory for every Gaussian in the grid: a 100 x 100 grid requires 10^ 
copies of the photometric data, yet it is already feasible if only the 
sectors are used. 

The advantage of this method is that it is guaranteed to find 
the lowest solution (although it may not be unique) within the 
accuracy imposed by the adopted grid of Gaussians. A big disad- 
vantage is that in practice many more Gaussians are used by the 
NNLS algorithm to produce a fit with the same as found with 
our nonlinear algorithm of Section ; 



3.3 



This is in general an impor- 
tant point, since one often needs to use the MGE parametrization to 
perform subsequent calculations whose speed depends on the num- 
ber of Gaussians: the computation of the potential scales linearly 
with the number of Gaussians, while for example the solution of 
the Jeans equations scales as A'^^ (Emsellem et al. 1994a). 

As an example we have used a 100 x 100 grid of (a, q') values 
(10* Gaussians) to produce an MGE model for NGC 4342 using the 
same input photometry of Section 3.3. The NNLS algorithm gave a 



best solution with 38 positive Gaussians and an error only slightly 
smaller than that of the 13 Gaussians nonlinear solution of Fig. ^ 
(RMS error of 1.7 per cent, to be compared with 1.8 per cent for 
the nonlinear solution). 

Since the NNLS solution is ill-conditioned it may be possi- 



ble to extract a subsample from the 38 Gaussians that gives es- 
sentially the same Xbest of the best solution. In particular one 
would like to find the smallest set of Gaussians that still verifies 
< (1 + Xbcst' where e represents a chosen fractional toler- 
ance. Unfortunately this is a combinatorial minimization problem 
whose exact solution would require (today) a very large amount of 
time even with 38 Gaussians. We have thus adopted the following 
approximate method: 

(i) solve the NNLS problem obtained by fitting the whole (a, q) 
grid of Gaussians Gj to the photometric data and define Xbcst = 

(ii) construct a set S containing the positive Gaussians Gj from 
the previous solution, with j = 1, . . . , P; 

(iii) for j = 1, . . . , P solve a new NNLS problem, excluding 
Gj from S, and record the Xj obtained in the solution. Define 

Xj,min 



rmn{x?, ■ • ■ ,Xp}; 



(iv) ifx'n.in < + 



2 

Xbcst ' 



, where e is the maximum accepted 
fractional error in , then eliminate Gj from S, update the number 
P ^ P — 1, and go back to step (iii), otherwise return S. 

Running the above fully linear algorithm, we have been able 
to find (in a time 400 x longer) a set of 14 Gaussians that produced 
a x^, in the NGC 4342 fit, not larger than that of the 13 Gaussians 
nonlinear solution. In general we have found that this linear algo- 
rithm is not competitive, both in speed and num ber o f Gaussians 
with the one described in Section 



3.3 



Its extreme 



for a given x 

robustness however is very attractive and it can be very useful to 
test that a given solution cannot be significantly improved, or in 
cases where the final number of Gaussians is not an issue. 

The solution obtained by the fully linear MGE fitting method 
can be further improved by the nonlinear fitting method and the 
steps (ii)-(iv) outlined above can be iterated once more. We have 
implemented all this and found that it works very well, although the 
extra complication of combining the two methods is not required 
most of the times. 

3.5 MGE fits of galaxies with isophote twist 



In Sections 5.3-3.4 we have shown how to fit MGE models to 



galaxies having isophotes with four-fold symmetry (constant PA 
with radius). Here we briefly discuss how to produce an MGE fit to 
galaxies whose isophotes twist with radius. 



The extension of the nonlinear fitting method of Section 3.3 to 
galaxies with isophote PA that varies with radius comes very natu- 
ral from the previous discussion. It has been shown that, in a sin- 
gle component object like an elliptical, every Gaussian contributes 
to the surface brightness essentially only at an elliptical radius m' 
close to its a. Consider for simplicity the case of constant axial ra- 
tio q' . As before the CTj will have to be logarithmically spaced and 
all the q'j can be set equal to a representative galaxy axial ratio. A 
very good starting guess for PA ipj of the Gaussian Gj will then 
be given by the PA of the isophote with semimajor axis a — Oj. 
The PA at the various radii can be measured before fitting with 
standard photometry packages. As before for the constant-PA fit- 
ting problem, the starting guess with constant axial ratio provides 
a good approximation to the actual model fit with varying isophote 
ellipticity and the nonlinear optimisation method converges to the 
global minimum. 

Fig. |ll] presents the profiles of the A'^ = 11 Gaussians MGE 
best-fitting model for a WFPC2/F702W image of the elliptical (E3) 
galaxy NGC 5831, that presents a ~ 35° isophote twist in the range 
- 40" (Peletier et al. 1990; Binney & Merrifield 1998, figure 
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Figure 11. Same as in Fig. ^ for the = 11 Gaussians MGE fit to tiie 
WFPC2/F702W pliotometry of NGC 5831. Only 6 equally spaced repre- 
sentative profiles are shown, of the A^sec = 36 sectors actually used in the 
fit. 0° conesponds to the PA of the central isophotes. The RMS error is 1.2 
per cent. 



4.37). The profiles have been measured along 36 sectors 5° wide, 
covering the whole galaxy image from —90° to 90°. Since we are 
fitting a centre- symmetric model, the counts measured within sec- 
tors symmetric with respect to the centre have been averaged before 
fitting. A very good MGE model can be obtained for this galaxy. As 
expected, although the PA is different for the various fitted compo- 
nents, the Gaussians succession is still regular and predictable as 
in the case of constant PA of Fig. A comparison between the 
contours of the data and of the convolved MGE model is shown in 
Fig.0. 

Fig. O shows an MGE fit to the central region of the SO barred 
galaxy NGC 2950. This object presents a sharp ~ 60° isophote 
twist inside the central 20", due to the presence of a double bar 
component (Friedli et al. 1996; Rest et al. 2001). AnN = 7 Gaus- 
sians MGE model is able to reproduce very well the central surface 
brightness of this object, including the strongly twisted and non el- 
liptical isophotes, with the exception of the slightly "spiral-like" 
isophote shape that peaks at an elliptical radius of ~ 7" and is 
clearly visible from the contour plot. 

The fully-linear MGE fitting method of Section 



3.4 can also 




Figure 12. Contour map for the 160" X 160" (1000 s) WFPC2/F702W im- 
age of elliptical galaxy NGC 5831. Note the ^ 35° isophote twist of this 
galaxy. Overlaid are the contours of the MGE best-fitting model whose pro- 
files are shown in Fig. 



formed with a three-dimensional grid (a, q' , ^p) of Gaussians. The 
only problem of this approach is that the memory requirements will 
increase considerably, but a 30x30x30 grid of Gaussians can al- 
ready be solved with current averages PCs. 



4 AVAILABILITY 

The complete IDL^ source code implementing the linear and non- 
linear MGE fitting algorithms, described in this paper, both with 
and without isophote twist, together with manual and usage exam- 
ples, can be found at http://www.strw.leidenuniv.nl/~mcappell/idl. 



be easily extended to deal with isophote twisting. It is simply nec- 
essary to add Gaussians with a range of PA ip, linearly spaced from 
—90° to 90°, in the representative library of Gaussian, for every 
{a, q') combination. The NNLS fit to the photometry will be per- 



5 CONCLUSIONS 

The MGE method is one of the few simple parametrizations that 
are general enough to reproduce the surface brightness of realistic 
multicomponent galaxies. In addition many dynamical and photo- 
metric quantities can be evaluated easily and accurately when the 
density is expressed in MGE form. 

In this paper we have described a simple yet powerful algo- 
rithm that reduces the process of generating an MGE fit to multiple 
galaxy images to a simple, fast and automatic task. We have pro- 
vided examples of its practical use and have tested it by accurately 
reproducing the photometry of a relatively large sample of galax- 
ies, both with and without isophote twists. We have also compared 
some of our fits with previously obtained photometric models. 

We have also described an alternative algorithm that, although 
currently less practical, due to the larger computing power re- 
quirements, is guaranteed to converge to the minimum solution 



http://www.rsinc.com 
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Figure 13. Contour map of the central 34"x34" of the (730 s) 
WFPC2/F814W (7-band) image of the barred SO galaxy NGC 2950. The 
best-fitting MGE model is overlaid. Note the sharply twisted and non ellip- 
tical isophotes due to a double bar within this galaxy. N = 7 Gaussians 
have been used in the MGE model and the resulting RMS en'or is 2. 1 per 
cent. 

within the accuracy imposed by an adopted grid of parameters in 
the solution space. 

These algorithms have been implemented in an IDL program 
that can produce an MGE model starting from the observed images 
of a galaxy, requiring the user to only input the coordinates of the 
galaxy centre, the PA and a characteristic flattening. Multiple res- 
olution images (e.g., ground-based and HST) can easily be fitted 
together, in one single step. The complete IDL source code imple- 
menting the algorithms described in this paper is made publicly 
available. 
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